load motor.mat;
XtestNorm = motorNormalize(Xtest);
XtrainNorm = motorNormalize(Xtrain);

X = linspace(-1, 1);
X = X';

meanPredictions = zeros(6,2);
meanPredictions(1,1) = 5;
meanPredictions(2,1) = 10;
meanPredictions(3,1) = 15;
meanPredictions(4,1) = 20;
meanPredictions(5,1) = 25;
meanPredictions(6,1) = 30;

RMSEtest = zeros(6,2);
RMSEtest(1,1) = 5;
RMSEtest(2,1) = 10;
RMSEtest(3,1) = 15;
RMSEtest(4,1) = 20;
RMSEtest(5,1) = 25;
RMSEtest(6,1) = 30;

RMSEtrain = zeros(6,2);
RMSEtrain(1,1) = 5;
RMSEtrain(2,1) = 10;
RMSEtrain(3,1) = 15;
RMSEtrain(4,1) = 20;
RMSEtrain(5,1) = 25;
RMSEtrain(6,1) = 30;

meanPredictions = zeros(26,2);
RMSEtest = zeros(26,2);
RMSEtrain = zeros(26,2);

W = zeros(31, 26);

for a = 1:26
   meanPredictions(a,1) = a+4;
   RMSEtest(a,1) = a+4;
   RMSEtrain(a,1) = a+4;
end


numDataPoints = size(XtrainNorm,1);
for Mindex=1:26
    mu = linspace(-1,1,meanPredictions(Mindex,1));
    sigma = 3*(mu(2)-mu(1));

    phi = ones(size(XtrainNorm,1),meanPredictions(Mindex,1)+1);
    for x=1:size(XtrainNorm,1)
        for j = 2:meanPredictions(Mindex,1)+1
            phi(x,j) = exp(-(XtrainNorm(x)-mu(j-1))^2 / 2*sigma^2);
        end
    end
    
    powerIndex = meanPredictions(Mindex,1);
    designMatrix = ones(numDataPoints,powerIndex+1);
    for rowDesign=1:numDataPoints
        for colDesign=2:powerIndex+1
            designMatrix(rowDesign,colDesign)=exp(-(XtrainNorm(rowDesign,1)-mu(colDesign-1))^2 / 2*sigma^2);
        end
    end
    w = (designMatrix' * designMatrix)\designMatrix' * Ytrain;
    
    YPredictTrain = zeros(size(XtrainNorm,1),1);
    for xtestIndex=1:size(XtrainNorm,1)
        YPredictTrain(xtestIndex,1) = phi(xtestIndex,:)*w;
    end
    
    
    testPhi = ones(size(XtestNorm,1),meanPredictions(Mindex,1)+1);
    for x=1:size(XtestNorm,1)
        for j = 2:meanPredictions(Mindex,1)+1
            testPhi(x,j) = exp(-(XtestNorm(x)-mu(j-1))^2 / 2*sigma^2);
        end
    end
    
    YPredict = zeros(size(XtestNorm,1),1);
    for xtestIndex=1:size(XtestNorm,1)
        YPredict(xtestIndex,1) = testPhi(xtestIndex,:)*w;
    end
    meanPredictions(Mindex,2) = sum(YPredict)/size(XtestNorm,1);
    
    %part 3 repeat
    for i=1:size(Ytest,1)
        RMSEtest(Mindex,2) = RMSEtest(Mindex,2) + (Ytest(i,1) - YPredict(i,1)).^2;
    end
    RMSEtest(Mindex,2) = sqrt(RMSEtest(Mindex,2)/size(Ytest,1));
    
    for i=1:size(Ytrain,1)
        RMSEtrain(Mindex,2) = RMSEtrain(Mindex,2) + (Ytrain(i,1) - YPredictTrain(i,1)).^2;
    end
    RMSEtrain(Mindex,2) = sqrt(RMSEtrain(Mindex,2)/size(Ytrain,1));
    
    w(31,1) = 0;
    W(:, Mindex) = w;
end
%figure(251);
%plot(meanPredictions(:,1),meanPredictions(:,2));

figure(252);
subplot(1,2,1)
plot(RMSEtest(:,1),RMSEtest(:,2));
subplot(1,2,2)
plot(RMSEtrain(:,1),RMSEtrain(:,2));

